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We present a local order parameter based on the standard Steinhardt-Ten Wolde approach that is capable both of tracking 
and of driving homogeneous ice nucleation in simulations of all- atom models of water. We demonstrate that it is capable 
of forcing the growth of ice nuclei in supercooled liquid water simulated using the TIP4P/2005 model using overbiassed 
umbrella sampling Monte Carlo simulations. However, even with such an order parameter, the dynamics of ice growth in 
deeply supercooled liquid water in all-atom models of water are shown to be very slow, and so the computation of free 
energy landscapes and nucleation rates remains extremely challenging. 
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I. INTRODUCTION 

It is well-known that substances cooled below their thermo- 
dynamic freezing point do not necessarily freeze, especially 
when they are very pure. Homogeneous nucleation is a kinet- 
ically disfavoured process; according to classical nucleation 
theory, a critical cluster must spontaneously form in the super- 
cooled liquid before crystallisation can proceed. ^"^ The kinetic 
barrier, in the framework of classical nucleation theory, arises 
from a competition between a favourable bulk free energy dif- 
ference between the phases and the unfavourable formation of 
an interface between the phases. The presence of a free energy 
barrier to nucleation makes homogeneous nucleation a rare 
event. 

Homogeneous nucleation has been studied using computer 
simulations in a range of systems however, the nucleation 
of ice, despite being of fundamental interest,-^ is still prov- 
ing to be difficult to simulate successfully despite the appar- 
ent simplicity of the process. The crystallisation of ice has 
been studied in a large number of simulations neverthe- 
less, whereas several simulations of homogeneous nucleation 
using the MW coarse-grained water model^^ have been rea- 
sonably successful,^^'^^'^^'^^ all-atom simulations have been 
less so. The MW potential is a good representation of the 
structure and the thermodynamics of water; however, it has un- 
realistically fast dynamics and no representation of hydrogens, 
and so the use of an all-atom model would offer considerable 
further insight into the process. However, Matsumoto and 
co-workers' single MD trajectory of TIP4P water nucleating 
into ice remains the only successful brute-force simulation of 
homogeneous ice nucleation with an all-atom model. ^ Other 
simulations have used small system sizes or looked at condi- 
tions that are not representative of homogeneous nucleation 
from the bulk liquid water. While brute-force simulations of a 
rare event are unlikely to be successful, the use of rare event 



methods can potentially allow us to compute free energy land- 
scapes for nucleation. Such calculations have recently been 
attempted in the homogeneous ice nucleation simulations of 
Radhakrishnan and Trout, who used umbrella sampling, ^^'^^ 
and Quigley and Rodger, who used metadynamics.^^ However, 
the use of global order parameters in driving homogeneous ice 
nucleation can lead to non-physical nucleation pathways, as 
we discuss below. 

It is a surprising state of affairs that modern simulation meth- 
ods have so far not been able to capture the fundamental phys- 
ical behaviour of the homogeneous nucleation process of ice. 
Some of the outstanding problems are how to simulate the ho- 
mogeneous nucleation of ice using a local measure of order, 
and the determination of a free energy landscape using such 
an order parameter. Here, we address one of these aspects: 
we develop some appropriate order parameters to allow us to 
drive the nucleation process and grow a single ice cluster using 
all-atom models of water. We first examine order parameters 
used in nucleation studies in general (Section II) and then in- 
troduce the order parameters we use in all-atom water model 
simulations (Section III). We discuss the water potential and 
simulations methods we used in Section IV, and we present 
the results of driving nucleation in Section V. Finally, we dis- 
cuss the outstanding problems that need to be overcome in 
order to obtain a free energy landscape and nucleation rates in 
Section VI. 

II. ORDER PARAMETERS IN NUCLEATION STUDIES 

In order to monitor the process of nucleation, we require 
a quantitative measure that can distinguish how far along the 
process is. The quantity describing this is usually known as 
an order parameter. The first step in deciding on an order 
parameter is to classify particles as being solid-like or liquid- 
like in nature. In the nucleation literature, such classification is 
often based on the Steinhardt classification parameter^^'^^ 
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Ylyn {Oij^ (pij) are the spherical harmonics, 6 and (p are the polar 
angles measured in an arbitrary laboratory frame of reference 
and A^neighs(0 is number of neighbours of particle i^^ It is 
important to note that all qi{i) are rotationally invariant regard- 
less of the choice of /. There is no radial component in this 
scheme; one can be introduced if necessary. However, a limited 
radial dependence arises through our definition of neighbours. 
In nucleation studies, we are often hoping to compare the envi- 
ronment about a particle to a symmetrical crystalline system, 
and so we can simply pick the value of / that best corresponds to 
the symmetry of the crystalline system and compute the spheri- 
cal harmonic expansion coefficients for only that /.^^ It is both 
convenient and computationally less expensive to replace the 
complex spherical harmonics with their real analogues 
we continue to denote complex conjugates in the following for 
generality, but they can be dropped if real spherical harmonics 
are used. 

A. Global order parameters 

In many studies, the local Steinhardt classification param- 
eters (Eq. (2)) are averaged across the system to give global 
Steinhardt order parameters; these are typically expressed as 
the magnitude of the vector sum of the local classification pa- 
rameters averaged over all N particles in the system, namely^^ 
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Although local classification parameters qi are a good measure 
of the local order round a particle, they are generally non-zero 
both in the solid phase and in the liquid phase, as the liquid 
is often reasonably well-ordered, especially when considering 
only its first neighbour shell. However, the vectors add inco- 
herently in the liquid phase and the global order parameter Qi 
averages out to zero for large systems, whilst it does not do so 
in the solid phase.^^ As a result, the increase in the magnitude 
of Qi can be used to track how solid-like a system is. Global 
order parameters of this type have been used in two previous 
studies of homogeneous ice nucleation. ^^'^^ However, global 
order parameters are not ideal in nucleation studies, particu- 
larly in studies where the system is not only tracked, but driven 
to increase the value of a global order parameter.^^ 

First of all, in nucleation studies, we typically wish to per- 
form calculations in such a way as to enable us to compare the 
results to classical nucleation theory. To do this, we need to 
know the size of the largest crystalline cluster, usually by know- 
ing how many particles there are in the cluster. However, with a 
global order parameter, we have no knowledge of what the size 
of the largest crystalline cluster is; indeed, the physical mean- 
ing of any particular value of a global order parameter is not 
only system-size dependent, but physically hard to interpret, 
and a free energy landscape calculated as a function of a global 



order parameter does not have a clear physical interpretation 
in terms of nucleation. 

Secondly, the surface free energy is in competition with 
the more favourable entropy arising from a larger number of 
smaller clusters. It can be shown^^ that this entropy can play 
a significant role in small systems that can be simulated on 
computers: in the early stages of nucleation, many small nuclei 
are always more favourable than one single nucleus comprising 
the same number of particles. There is, however, a crossover 
to the expected behaviour once a certain nucleus size has been 
passed. This means that, when comparing results to classical 
nucleation theory, a 'global' measure of crystallinity - which 
effectively induces an entropic break-up of small clusters - is 
not sufficient.^^ 

Finally, we have suggested in our previous work^^ that the 
pathways produced when the nucleation process is driven by 
global order parameters may be inconsistent with the natural 
nucleation pathways, particularly so in the case of ice. Since 
there is no distinction between the liquid and the solid states 
of particles when global order parameters are used, driving the 
system to increase its global order parameter can induce an 
orientational coherence even in the liquid state. This suggests 
that a particle is influenced not only by its neighbours, but 
potentially by a crystalline cluster that is very far removed from 
it, even though such a long-range interaction has no basis in 
reality.^^ Furthermore, when rare event techniques are applied 
to a system, it is relatively easy to compensate for arbitrarily 
large free energy barriers; there is a danger, therefore, that if 
the natural, lowest free energy pathway is dynamically slow 
(having accounted for the free energy barrier associated with 
the process itself), such a pathway may not be observed and a 
higher free energy pathway could be found instead provided 
that it is dynamically faster and its higher free energy has 
been negated by a rare event method. We have previously 
also suggested that, for ice nucleation, the use of global order 
parameters in driving nucleation may lead to precisely such 
non-physical high free energy pathway s.^^ 

B. Local order parameters 

Even though the local Steinhardt classification parameters 
defined above cannot distinguish between solid and liquid par- 
ticles on their own, there are a few approaches that allow us 
to do this without sacrificing their local nature. It is often 
convenient to calculate the dot products of the individual lo- 
cal classification parameters expressed in vector form, 
whose (2/ + 1) components are the Steinhardt parameters qim{i) 
for m G [—/,/] n Z, with the equivalent vectors of a particle's 
neighbours. We then calculate the rotationally invariant func- 
tion di{ij) = Qi (/) • qfij), where / and j are neighbours; this 
dot product value ranges between —1 and +1. By plotting 
the distribution of J/ values for the liquid and the crystalline 
phases, a critical threshold dc can be determined as the first 
point where the probability of being in the solid phase is non- 
zero. The number of crystalline connections is then defined 
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where H is the Heaviside step function. The number of con- 
nections should be higher in the soHd phase than in the hquid 
phase, and a criterion involving a threshold minimum number 
of connections to distinguish between the two phases is often 
a good classification parameter. ^^'^^ 

Another procedure involves calculating the neighbour- 
averaged contribution,^^'^^ 
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where j runs over all the neighbours of particle /, and includes 
the particle itself (when j = 0). The average local bond classifi- 
cation parameter is then given by 
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In both approaches mentioned, the second neighbour shell is 
effectively taken into account through the use of local Stein- 
hardt vectors of the first neighbour shell, either by averaging 
or by taking dot products. 

Any two particles belong to the same crystalline cluster if 
they are both classified as being solid-like and are located 
within a certain fixed distance of one another (that is, they 
are neighbours). Once all the particles have been classified, 
the size of the largest such cluster is normally calculated in 
nucleation studies; this then acts as the overall (local) order 
parameter. 

It is very useful for an order parameter used in driving nucle- 
ation to be local in nature; however, that an order parameter is 
local is not sufficient for it to be a valid metric used for driving 
the nucleation process. For example, in their ice nucleation 
simulations,^ ^'^^ Brukhno and co-workers used a maximum 
director projection approach to yield local order parameters. 
However, although these order parameters do permit the growth 
of ice to be driven in a fixed orientation with respect to the 
simulation box, the rotational bias inherent in the procedure 
induces a non-local and non-physical orientational coherence 
in the growing ice cluster; we suggest that, as a result, this 
order parameter is not suitable to study ice nucleation. 

More complex order parameters used to track larger 
molecule crystallisation have also been proposed,^^ whilst an 
overview of many simpler order parameters used in vapour- 
liquid nucleation was produced by Senger and co-workers.^ ^ 

III. ORDER PARAMETERS FOR HOMOGENEOUS ICE 
NUCLEATION 

The choice of an order parameter to drive ice nucleation is 
not trivial: as discussed above, it is preferable that it be local; 
it must be forgiving enough to be able to induce the growth 
of a small ice cluster; and it must be strict enough to ensure 
that the structure grown is actually ice-like and has, ultimately, 
long-range order. In our simulations, we use a variation of 
the dot product approach described above that is commonly 
used in studies of tetrahedral liquids. ^^'^^'-^^'"^^ We choose to 
use 1 = 3, since the / = 3 spherical harmonics are the ones best 
describing tetrahedrality. A plot of the distribution of d^{ij) is 
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FIG. 1. A typical probability density distribution for all pairs of 
^3 (iJ) = Q3 (0 ■ ^3 (/)' where the centres of mass of molecules / and 
j are within 3.5 A of each other. The three states depicted were equi- 
librated at 200 K (using the TIP4P/2005 water model) and the ice 
structures are not, therefore, 'perfect'. This figure is analogous to 
those in Refs 25, 36, and 52. 

shown in Fig. 1; to account for the eclipsed bond in hexagonal 
ice, we define a classification parameter as 

A^neighs(0 
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where 



^ fl if [(x< -0.825) V(-0.23<x< 0.01)], 
1 otherwise. 



These limiting values were chosen to encompass d^iij) re- 
gions (Fig. 1) where the probability density function for either 
ice phase has a value greater than 0.1. We classify a molecule 
as ice-like if /tconnections ^ 3 and as liquid-like otherwise. This 
gives perfect identification in both equilibrated cubic and equi- 
librated hexagonal ice. The order parameter we use to track 
the progress of nucleation is the size of the largest cluster of 
molecules classified as ice, where two molecules belong to 
the same crystalline cluster if they are both ice-like and their 
centres of mass are within 3.5 A of each other. 

Unfortunately, whilst this classification procedure (regard- 
less of the precise details of the parameterisation of the limiting 
values) works well for the MW potential,^^'^^ it does not do 
so in all-atom models of water. When used in the form pre- 
sented above in umbrella sampling^^ or forward flux sampling 
simulations,^^ natural fluctuations in the system often result 
in molecules satisfying the order parameter even if they are 
not really ice-like. When forced to grow with a bias sing po- 
tential, 'chains' form more easily than real ice grows, even 
though such chains actually represent an abuse of the order 
parameter, and 'ice' structures as depicted in Fig. 2 are com- 
monplace. The real issue is not just that such chains form, but 
that when they do form, the system is not subsequently able to 
grow the correct ice structure, and ice growth is arrested. In 
order to alleviate this problem, we (a) classify any molecule 
with more than four neighbours (within 3.5 A) as being liq- 
uid, and (b) explicitly exclude molecules belonging to chains 
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FIG. 2. An example of non-physical chain growth in TIP4P/2005 
umbrella sampling simulations when using the order parameter with- 
out chain removal as described in the text. The system has 1000 
molecules at 240 K, starting from a 24-molecule cluster. Molecules 
classified as being part of the largest crystalline cluster are shown in 
red and violet; there are 45 molecules in this cluster. Molecules whose 
centres of mass are within 3.5 A are connected with lines. Molecules 
shown in violet would be removed from the largest crystalline cluster 
on application of the chain removal algorithm described in the text. 

from the largest crystalline cluster. We achieve the latter by 
removing any molecules with only one neighbour belonging to 
the largest cluster from the largest cluster, except if that single 
neighbour is connected to three further molecules in the largest 
cluster. This allows 'chains' comprising a single molecule to 
form and thus allows ice to grow. We iterate the procedure 
until no further molecule is removed. 

Removing chains could be problematic in the initial stages 
of nucleation: it is impossible for an ice structure smaller than 
a single chair (or boat) not to be formed of chains, and small 
rings may form instead if forced. However, due to the similar- 
ity of liquid water and ice, it is possible to wait for a boat or 
chair cluster to form spontaneously in a simulation, and limit 
umbrella sampling to systems with clusters larger than ~10 
molecules. Doing so does not appear to be a significant limita- 
tion when the critical cluster is expected to comprise over 100 
molecules. 

Although neighbour-averaged classification parameters are 
excellent at distinguishing between the phases of systems with- 
out much structure in the liquid phase,^^'^^ it is less clear 
whether the same applies to well- structured liquids like wa- 
ter. Using / = 4 and / = 6 (as depicted in Fig. 3) gives better 
separation between the phases than does using 1 = 3. Liquid 
water and hexagonal ice are less well- separated in the {q4)-{qe) 
plane than in the Lennard- Jones case; nonetheless, a choice 
of {qe{i)) > 0-7 as an ice-liquid boundary would appear to 
be reasonable. Calculating the size of the largest cluster with 
this method results on average in only slightly smaller clusters 
compared to those resulting from the dot product approach, 
and chains are as problematic as when using the latter. We 
have shown^^ that for the MW monatomic model of water,^^ 
a dot product approach leads to a free energy landscape that 
is almost entirely consistent with classical nucleation theory. 
By contrast, including 'surface' molecules in the largest ice 
cluster, as attempted in other studies with the same model 
of water,^^'^^ appears to reduce the agreement with classical 
nucleation theory. Lechner and co-workers have recently 
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FIG. 3. Neighbour- averaged order parameters for systems of ices Ih 
and Ic and liquid water. All systems were equilibrated at 200 K using 
the TIP4P/2005 water model, and they contain different numbers of 
molecules. The neighbour cutoff distance was 3.5 A. 

used a combination of dot product vectors and their neighbour- 
averaged classification parameters to study the role of the sur- 
face and the bulk terms when comparing simulation data to 
classical nucleation theory for a soft-core colloid model.^^'^^ 
Using neighbour-averaged classification parameters reduces 
the size of the critical cluster when compared to one calculated 
using a dot product approach in their work.-^^'-^^ The clusters 
we have analysed with both approaches show only a small 
difference in cluster size; on average, the neighbour- averaged 
clusters are slightly smaller, but for individual configurations, 
the converse can also hold. Nonetheless, given that the number 
of molecules classified as being part of the cluster by the dot 
product approach is itself a rather conservative estimate of the 
cluster size,^^ and since it would be computationally extremely 
expensive to perform two-dimensional umbrella sampling with 
an additional order parameter, we restrain ourselves, for the 
time being, to using the dot product approach only. However, 
the neighbour-averaged approach is an attractive alternative, 
and ensuring that clusters are of essentially the same size with 
both approaches is a useful confirmation that the exact details 
of the order parameter do not appear to change the outcome 
significantly. 

IV. SIMULATION DETAILS 

An empirical model that reproduces a large number of ex- 
perimental results at a reasonable computational cost is the 
TIP4P/2005 model,^^ which seems to be the best of the 'sim- 
ple' all-atom models available at present, and is one that works 
well across several phases^^"^^ and at large supercoolings.^^ 
Although more complex models can, at a considerable compu- 
tational expense, capture more of the underlying physical be- 
haviour, their use, at least in their present state of development, 
does not necessarily result in a better description of water. ^^"^^ 
In this work, we therefore use the TIP4P/2005 model of water, 
as described in the original paper by Abascal and co- workers, 
using the same parameters for cutoffs and Ewald summation. 

To simulate the nucleation of ice, we use the Metropolis 
Monte Carlo (MC) approach^^ in the isobaric-isothermal en- 
semble, coupled with umbrella sampling^^ to drive the pro- 
cess. Our umbrella sampling simulations use adaptive um- 






FIG. 4. Representative nucleation snapshots from umbrella sampling simulations of TIP4P/2005 water. In each case, two pictures depict the 
same cluster from different perspectives; one within the liquid framework (in cyan) and one showing solely the largest crystalline cluster. In 
the former, spheres represent centres of mass of molecules classified as ice: red spheres correspond to cubic ice, orange spheres correspond to 
hexagonal ice and pink spheres correspond to ice molecules not within the largest crystalline cluster. Pictures representing solely the largest 
cluster depict both the oxygen (red) and the hydrogens (white) of each molecule. In (a), an 82-molecule ice cluster grown from the supercooled 
liquid at 240 K is shown; in (b), a 73-molecule ice cluster grown from a small cluster of Ih ice at 240 K is shown; and in (c), a series of ice 
clusters of increasing size (comprising 23, 60, 77, 107 and 145 molecules from left to right) grown from a small cluster of Ic ice at 200 K is 
depicted. There are 1900 molecules in the system in (a) and the first three configurations of (c), and 2500 molecules in (b) and the last two 
configurations of (c). Simulations of nucleation from a hexagonal seed (shown in (b)) were undertaken using the hybrid Monte Carlo approach, 
and the rest by a standard Metropolis Monte Carlo approach, p = 1 bar. 



brella weights^^ rather than the more traditional quadratic bias 
potentials.^ ^ 

All-atom models of water have particularly slow dynamics, 
especially for crystal growth, making any possible improve- 
ment in computational speed worth considering. It has been 
shown that ice crystal growth occurs more rapidly (in computer 
time) in molecular dynamics (MD) simulations than in corre- 
sponding Monte Carlo ones,^^'^^ which may suggest that the 
collective motion possible in MD simulations helps to speed 
up the dynamics of cluster reorganisation, and thus aids us in 
driving crystallisation. The choice of MD simulations over 
MC simulations is therefore appealing; however, as we are ulti- 
mately interested in the free energy landscape of ice nucleation, 
we wish to continue to use umbrella sampling and thus Monte 
Carlo simulations. 

While the majority of simulations presented in this paper 
used the standard Monte Carlo method, we have recently be- 
gun to couple Monte Carlo with MD simulations in a hybrid 
Monte Carlo approach,^^ where short MD simulations replace 
single particle Monte Carlo moves. Provided that the MD in- 
tegrator is time reversible and symplectic and that the choice 
of momenta from the Maxwell-Boltzmann distribution is ac- 
counted for in the Metropolis acceptance criterion, detailed 
balance is obeyed^^ irrespective of the fact that the MD hamil- 



tonian does not incorporate an umbrella sampling term. We 
implement the symplectic and time reversible quaternion-based 
algorithm of Miller III and co-workers^^'^^ to simulate rigid 
body rotations. We are able to drive nucleation using both 
methods; however, simulations are considerably faster in real 
time when the hybrid Monte Carlo approach is used, which 
confirms the importance of collective motion for nucleation. 

V. NUCLEATION PATHWAYS 

We have run umbrella sampling simulations in a variety of 
systems with TIP4P/2005 water. Typical simulations involved 
between 1900 and 2500 water molecules and three distinct 
scenarios were considered: growth from a seed hexagonal ice 
cluster, growth from a seed cubic ice cluster and growth directly 
from the supercooled liquid water. 

The approach we used was to bias the umbrella sampling 
weights to favour larger clusters, although the weights cho- 
sen were only slightly overbiassed compared to the classical 
nucleation theory prediction of the free energy barrier. If the 
umbrella weights bias the growth to be too quick, then the 
system can begin to grow defective crystal nuclei that cannot 
repair themselves by shrinking and regrowing: the weights 
must be sufficiently small to allow clusters to grow and shrink 
throughout each umbrella sampling window. 
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Several snapshots of ice growth in such simulations are de- 
picted in Fig. 4. Provided that umbrella sampling does not 
attempt to drive the nucleation too quickly, the resulting ice 
clusters appear to be very reasonable. This suggests that the 
order parameter presented above is a suitable order parame- 
ter both to track and to drive the process of homogeneous ice 
nucleation. 

It is intriguing to note that it appears that in what we believe 
to be overbiassed driving of nucleation, cubic ice seeds seem 
to grow in a cubic fashion, and conversely, hexagonal ice seeds 
result in the growth of hexagonal ice. Growth directly from 
the supercooled liquid is somewhat more tricky: because the 
order parameter we use cannot track clusters smaller than 6 
molecules, and cannot readily track the growth of very small 
clusters of ice, it was necessary to wait for small clusters to 
form spontaneously. These were then biassed to grow fur- 
ther, although with a very gentle set of umbrella weights, and 
only when the growth became spontaneous with that set of 
weights did we progress to higher umbrella sampling windows. 
Analogously to what we observed in MW simulations,^^ ice 
grown directly from the supercooled liquid contains both cubic 
and hexagonal ice patterns, with cubic ones dominating, but 
less so than in the corresponding MW clusters: for example, 
60-molecule ice clusters in this work were classified to have 
approximately 70 % core cubic ice, whilst the MW analogues 
were about 90% cubic. Whether this is a result of a true dif- 
ference between the TIP4P/2005 and MW models of water or 
simply a consequence of overbiassed non-equilibrated driving 
in this work is unclear and warrants further investigation. The 
ice clusters observed in this work are roughly spherical, and 
this sphericity follows the same trends as for the MW model 
nucleation reported previously. 

Although we made no serious attempt to equilibrate the nu- 
cleation simulations, it may nevertheless be possible to gain 
physical insight into the nucleation pathway from our simu- 
lations. In order to compare the pathway observed with that 
reported in our simulations using the MW potentiaP^ and with 
Quigley and Rodger's simulation of TIP4P nucleation,^^ we 
calculate the Steinhardt- style and Chau-Hardwick-style^^ 
tetrahedrality parameters as defined by Quigley and Rodger,^^ 
including their smoothing function. These are given by 
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where is the number of particles and /(?■,)) is the smoothing 
function, and 
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where fij is the unit vector from particle / to particle j. The 
smoothing function is defined as 
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FIG. 5. The global order parameters ^^d ^ calculated as a function 
of the size of the largest crystalline cluster, the order parameter used 
to drive nucleation in this work, for the system seeded with a cubic 
ice nucleus. Error bars show the standard deviation for the population 
of configurations at each cluster size. The results depicted here refer 
to the 576 particles nearest the centre of mass of the ice nucleus. 
r = 200K,;7 = Ibar. 

We take into account the nearest 576 molecules from the centre 
of mass of the largest crystalline cluster as determined by the lo- 
cal order parameter defined above in order to ensure that these 
results can be compared to the previous work. The resulting 
diagram for the simulation of nucleation from a cubic ice seed 
nucleus is depicted in Fig. 5. Although the curves are rather 
noisy, as can be expected from a set of simulations that have 
not been equilibrated, and there are clear minor variations in 
slope corresponding to different umbrella sampling windows, 
we can nevertheless observe that the two order parameters plot- 
ted change roughly linearly as the cluster size increases. This 
linearity, which reflects the growth of an ice nucleus into a 
largely unperturbed liquid, is consistent with the nucleation 
pathway we reported for the MW model of water,^^ although 
the actual values of the global order parameters suggest that 
the system studied here is less well ordered than its MW ana- 
logue. This is perhaps not surprising considering that the MW 
systems we studied previously were very well equilibrated. Im- 
portantly, the pathway is rather different from that observed in 
Qe-l^ space by Quigley and Rodger,^^ as their free energy land- 
scape involves an initial increase in the orientational order 
before the tetrahedrality parameter begins to change, which 
implies that the entire system, rather than just a crystalline 
nucleus, becomes more ordered prior to the nucleation event. 
This adds further weight to our contention that global order 
parameters may locate pathways that are not fully consistent 
with the natural nucleation pathways, probably because such 
order parameters can induce global coherence in the system. 

VI. DISCUSSION AND CONCLUSIONS 

We have introduced an order parameter that is capable of 
tracking and driving the homogeneous nucleation of ice with 
the TIP4P/2005 water model. The order parameter is local 
in nature and thus does not exhibit the anomalous behaviour 
associated with global order parameters. We believe it to be the 
first rotationally invariant order parameter of this kind that is 
capable of driving homogeneous ice nucleation in simulations 
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of an all-atom model of water. One of the major difficulties 
in the development of such an order parameter is that the time 
that is required to confirm whether an order parameter is fit for 
purpose is very significant, given that the dynamics at reason- 
able supercoolings are so slow. In particular, it is important that 
the umbrella sampling weights not be increased too quickly 
even if it appears that no growth is forthcoming: one must 
exercise a considerable degree of patience when performing 
such nucleation simulations. 

The order parameters we have introduced represent the first 
step in obtaining free energy landscapes and nucleation rates 
for the homogeneous nucleation of ice from simulations. How- 
ever, there remain considerable challenges ahead. The basic 
problem is that the dynamics of ice nucleation for all-atom 
models are excruciatingly slow,^^'^^ making equilibration very 
difficult. This problem is illustrated by Fig. 6, which shows that 
at a reasonable 20% supercooling (200 K) in the TIP4P/2005 
model, a crystalline cluster of approximately 220 molecules in 
a system of 2500 molecules neither shrinks nor grows. Al- 
though only the first 1 ns of this pure MD simulation is shown 
in Fig. 6, the results are no different when simulated up to 70 ns, 
taking approximately a month of CPU time for each trajectory. 
Furthermore, to allow us to equilibrate a system, we would 
require not one, but a large number of freezing and melting 
events to be able to be simulated in such a period of computer 
time. 

One control parameter at our disposal is, of course, the tem- 
perature. The fastest rate of growth of an ice- water interface 
is generally some 10 K below the freezing point,^^"^^'^^ where 
the increased freezing driving force of cooler systems optimally 
balances their slower dynamics. The faster dynamics are illus- 
trated by the reasonably rapid melting of the crystalline cluster 
at 240 K (~5 % supercooling) as depicted in Fig. 6. One con- 
sideration that must be borne in mind when choosing a suitable 
temperature at which to perform simulations is that the majority 
of experimental rates have been reported at temperatures corre- 
sponding to supercoolings of between 27% and 10%;^^ these 
also encompass the atmospherically relevant conditions. While 
higher temperature simulations may be easier to run in some 
ways due to the expedited dynamics, raising the temperature 
is not without its problems. For example, while we can negate 
the higher nucleation free energy barrier with umbrella sam- 
pling, it is not just the barrier height, but also the critical cluster 
size that increases with increasing temperature. Increasing the 
temperature would thus require us to simulate considerably 
larger systems than are computationally affordable in order to 
avoid spurious finite size effects: for illustration, at 200 K, the 
critical cluster for TIP4P/2005 water is predicted by classical 
nucleation theory to encompass 10^ molecules, while at 240 K, 
this rises to 10"^ molecules. We could attempt to extrapolate 
the free energy barrier to lower temperatures using histogram 
reweighting^^ based on the results of small cluster simulations 
at higher temperatures. However, the calculation of nucleation 
rates requires the simulation of critical clusters, and so must be 
performed at sufficiently low temperatures so that the critical 
cluster is small enough to be feasible to simulate. 

How does one, then, successfully simulate the homogeneous 
nucleation process and obtain a free energy landscape and nu- 
cleation rate? The simplest strategy is simply to wait for a 
very long time: however, given how computationally challeng- 
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Time / ps 

FIG. 6. MD simulations of melting. The starting point is a crys- 
talline cluster comprising approximately 220 molecules embedded in 
supercooled liquid water. The curves exhibiting melting were sim- 
ulated at 240 K, whilst the remaining ones were simulated at 200 K. 
2500 TIP4P/2005 water molecules. Note that the melting point of 
TIP4P/2005 ice is 252K.61 

ing the process is, this may involve an inordinate amount of 
computer time. A second approach is to use more efficient 
simulation algorithms; indeed, as discussed above, the use of 
the hybrid Monte Carlo approach gives a significant advantage 
over standard Monte Carlo simulations in terms of simulation 
speed. Other tricks of the trade that might be advantageous 
include the use of hamiltonian exchange^^ to couple the sys- 
tem of interest to one that is dynamically faster, or the use 
of reaction fields in place of the computationally expensive 
Ewald summation. Finally, the water potential we use is an- 
other parameter of the system that is under our control. It 
has been suggested that there are few differences in the dy- 
namics of melting of most common all-atom ice models, 
and so a possible solution to the equilibration problem may 
be to use a model of water that is not necessarily the best at 
describing most experimental properties, but one whose dy- 
namics are computationally faster; examples might include 
TIP5P(-E)^^'^^ and the Nada-Van der Eerden potential. ^^'^^ 
Even though such potentials may be computationally more 
demanding than TIP4P- analogues on a per-step basis, the dy- 
namics of ice growth might nevertheless be faster. 

In conclusion, the development of a seemingly rigorous or- 
der parameter may help us to advance our understanding of 
ice nucleation; in particular, we have further corroborated our 
hypothesis that the difficulty in simulating ice nucleation in 
all-atom models such as TIP4P/2005 is more a result of the 
slow dynamics of the process rather than of an overwhelmingly 
large free energy barrier. We hope that the order parameter 
we have presented here represents a stepping stone towards 
the successful determination of a free energy landscape and 
nucleation rate for homogeneous ice nucleation for all-atom 
water models. 
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